title: “broccoli-population-structure”

author: “Sean O’Neill”

date: “8/2/2021”

output: html_document

html_document:

highlight: pygments

theme: united

toc: yes

toc_float:

collapsed: yes

smooth_scroll: yes

pdf_document:

toc: yes

rm(list=ls())

#get all the names
taxalabels <- read.delim("taxa/taxaBNUM_NAME.txt", header=F, stringsAsFactors=F, row.names = NULL)$V1
#make duplicate names unique
taxalabels <- make.names(taxalabels, unique = T)

#make some colors.   http://vrl.cs.brown.edu/color
clist <- list(
  sean=c("#fd00fe", "#62ce75", "#473c85", "#24a5f7", "#9a1073", "#4ad9e1", "#284e37", "#cad7d4", "#723521", "#fea53b", "#c00018", "#f7b5be", "#788c3b", "#5920af", "#eb74c3", "#6975fe", "#238910", "#528e8c", "#65f112", "#a57784", "#b3e61c", "#ae56ff", "#d6c951", "#ff7074"),
  shawn=c("#fd00fe", "#84bc04", "#1f3e9e", "#e775cc", "#46ebdc", "#256b33", "#82d1f4", "#7c28a4", "#2af385", "#214a65", "#aaa4e1", "#a1085c", "#148fae", "#cadba5", "#943112", "#fa8d80", "#604020", "#f1d438", "#eb1138", "#f7931e", "#270fe2", "#a07d62")
) 
# add length of palettes
lengths <- sapply(clist,length)
names(clist) <- paste0(names(clist),"_",lengths)

# par(mar=c(0.2,6,0.2,0))
# par(mfrow=c(length(clist),1))
# 
# #for some reason, barplot() doesn't work in .Rmd
# for(i in 1:length(clist))
# {
#   {barplot(rep(1,max(lengths)),col=c(clist[[i]],rep("white",max(lengths)-length(clist[[i]]))),axes=F,border=F)
#   text(x=-0.1,y=0.5,adj=1,label=names(clist)[i],xpd=T,cex=1.2)}
# }

Phylogenetic Tree

#source('scripts/rename-taxa.sh --infile="BOL21afZAXA00/RAxML/outraxml.raxml.supportFBP" --outfile="BOL21afZAXA00/RAxML/raxml.tree" --repfile="replaceNEWtoOLD.csv"')
#setwd("BOL21afZAXA00/")

out.group <-   c("W001_W001_WT_A","W006_W006_WT","W005_W005_WT","W004_W004_WT","W002_W002_WT","W003_W003_WT")

tree <- read.tree(file="RAxML/raxml.tree") -> uniqtree
tree$tip.label <- make.names(tree$tip.label, unique = T)
tree <- root(tree, outgroup = out.group)
tree <- ladderize(tree)


temp <-plot(tree,node.pos=1, font=1, align.tip.label=TRUE ,plot=FALSE)

is_tip <- tree$edge[,2] <= length(tree$tip.label)
ordered_tips <- tree$edge[is_tip, 2]
write.csv(tree$tip.label[rev(ordered_tips)], file="order.csv", row.names=FALSE)

# # Tree for popstructure
# svg("image/raxml.tree.svg", height = 35, width=20, pointsize = 10)
# par(oma=c(0,0,0,10))
# plot(tree,node.pos=1, font=1, align.tip.label=TRUE, cex=1.5, edge.width=4, no.margin=TRUE,  label.offset = 0.01 )
# dev.off()

tree.file <- "../raxml.tree_delete.png"
png(tree.file, units = "in",  height = 25, width=14, pointsize = 10, res = 150)

# tree.file <- "images/raxml.tree_delete.svg"
# svg(tree.file, height = 25, width=14, pointsize = 10)

par(oma=c(0,0,0,10))
plot(tree,node.pos=1, font=1,show.tip.label = TRUE,  align.tip.label=TRUE, cex=1, edge.width=2, no.margin=TRUE,  label.offset = 0.01 )
tree.labels <- tree$node.label
tree.labels[as.integer(tree.labels) > 0] <- "" 
#node.indices <- which(as.integer(tree.labels)<50)
#tree.labels <- as.character(tree.labels)
nodelabels(tree.labels, frame = "none", adj=c(-.125),cex=.75)
dev.off()  

knitr::include_graphics(tree.file)

PCA

#First, manually create a file called order.groups.csv. BOL2137SMinCn01 has an example

#choose a color list from clist.
gclist <- clist$sean_24

pca <- read.delim("PCA/PCA1.txt",skip = 2)
# pca2 <-read.delim("PCA/PCA2.txt",skip = 2)
# pca3 <-read.delim("PCA/PCA3.txt",skip = 0)
# hist(pca3$Eigenvector4)

pca$Taxa <- taxalabels
order <- read.delim("order.groups.csv")
#rearrange pca rows to be like order
pca <- pca[match(order$taxa,pca$Taxa),]
pca$group <- order$n00
pca$color <- gclist[order$n00]

#now make the group color list
gc <- order$color


axis = list(showline=FALSE,
            zeroline=FALSE,
            gridcolor=rgb(1,1,1), #'#ffff',
            ticklen=0,
            automargin = TRUE,
            titlefont=list(size=13))

fig11 <- pca %>%
  plot_ly(    
    height=1440,
    width=1440
    )
fig11 <- fig11 %>%
  add_trace(
    type = 'splom',
    
    dimensions = list(
      list(label='PC1', values=~PC1),
      list(label='PC2', values=~PC2),
      list(label='PC3', values=~PC3),
      list(label='PC4', values=~PC4),
      list(label='PC5', values=~PC5)

    ),
    text=~Taxa,
    diagonal=list(visible=F),
    marker = list(
      color = ~color,
      #colors = clist,
      size = 5,
      line = list(
        width = 1,
        color = 'rgb(1,1,1)'
      )
    )
  )
fig11 <- fig11 %>%
  layout(
    title = "Scatterplot Matrix for PCA",
    hovermode='closest',
    autosize = F,
    dragmode = 'select',
    plot_bgcolor='rgba(1,1,1, 0)',
    xaxis=axis,#list(domain=NULL, showline=F, zeroline=F, gridcolor='#ffff', ticklen=2, titlefont=list(size=13)),
    yaxis=axis,#list(domain=NULL, showline=F, zeroline=F, gridcolor='#ffff', ticklen=2, titlefont=list(size=13)),
    xaxis2=axis,
    xaxis3=axis,
    xaxis4=axis,
    xaxis5=axis,
    yaxis2=axis,
    yaxis3=axis,
    yaxis4=axis,
    yaxis5=axis

  )
fig11 <-  fig11 %>% style(showupperhalf = F)
fig11

Plot fastStructure

#setwd("BOL21afZAXA00/")


#Make a list of paths and a list of .meanQ file names
ffiles <- list.files("fastStructure",full.names=T, pattern = "*.meanQ")
flist <- readQ(files=ffiles)

#imports the order of taxa saved in the phylo tree steps.
tree.order.names <- read.delim("order.csv",header = T,)$x
#sorts the taxa the same way for pop structure plots.
tree.order.inds <- match(tree.order.names,taxalabels)


flist <- sortQ(flist, "k")
run.names <- attr(flist,which = "names")
if(length(unique(sapply(flist,nrow)))==1) flist <- lapply(flist,"rownames<-",taxalabels)
indfun <- function(x) flist[[x]][tree.order.names,]
flist <- lapply(run.names,indfun)
attr(flist, which = "names") <- run.names

#lab1 <- read.delim("TreeGroups000.txt", header = FALSE)

##Make the Plots!

text.size <- 5

kl <- length(flist)
kl1 <- round(kl/2)
kind1 <- c(1:kl1)
kind2 <- c((kl1+1):kl)

IchooseK <- c(8:12)

#uncomment this to pick your own values!
kind1 <- IchooseK

fplotfilename1 <- "fplot1"
fplotfilename2 <- "fplot2"

fastStructurePlot1 <- plotQ(alignK(flist[kind1], type = "across"), showlegend = T, imgoutput="join",splab = as.character(kind1+1), showindlab=T, useindlab=T, returnplot=T, exportplot=T,  basesize =  text.size, imgtype="png", exportpath = "images/", height=2,width=44, sppos="left",dpi=600,outputfilename=fplotfilename1)

## Drawing plot ...

## images//fplot1.png exported.

#fastStructurePlot2 <- plotQ(alignK(flist[kind2], type = "across"), showlegend = T, imgoutput="join",splab = as.character(kind2+1), showindlab=T, useindlab=T, returnplot=T, exportplot=T,  basesize = text.size, imgtype="png", exportpath = "images/", height=2,width=44, sppos="left",dpi=600,outputfilename=fplotfilename2)

#grid.arrange(fastStructurePlot1$plot[[1]])
#grid.arrange(fastStructurePlot2$plot[[1]])
#![Image Title](paste("images/",fplotfilename1,".pdf",sep="")){width=65%}

fm1 <- image_read(paste("images/",fplotfilename1,".png",sep=""))
image_rotate(fm1, 90) %>% image_write(paste("../../",fplotfilename1,"rot90",".png",sep=""))

#fm2 <- image_read(paste("images/",fplotfilename2,".png",sep=""))
#image_rotate(fm1, 90) %>% image_write(paste("../../",fplotfilename2,"rot90",".png",sep=""))

knitr::include_graphics(paste(fplotfilename1,"rot90",".png",sep=""))